Objective

Integrative analyses to assess the effect of baseline monocyte FACS, CSP IgG, and low annotation BTM expression on protection outcome by treatment (placebo or 1.8x10^6 PfSPZ). Model interactions.

Load libraries

library(tidyverse)
library(Biobase)
library(ggpubr)
library(pheatmap)
library(gtsummary)
library(flextable)
library(officer)
library(reshape2)
library(googledrive)
library(devtools)
source_url("https://raw.githubusercontent.com/TranLab/ModuleLists/main/Gene2ModuleExpressionScores.R")

Read-in data

Read-in data using google drive (may require authentication using gmail/google account)

#from google drive
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1NlFi_kCOrC2Hj2TkYgRMV3PCDELysSCI"), path = temp, overwrite = TRUE)
myDat_monocytes_CSPAb <- readRDS(file = dl$local_path)

Arrange data

basefont <- "sans"
basefontsize <- 10
pvalfontsize <- 4
mySignifLabel <- "p.signif"
symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("****", "***", "**", "*", "trend","ns"))


Mono_v_CSP_plot_nrow2_dat <- myDat_monocytes_CSPAb %>%
  dplyr::filter(treat == "1.8 x 10^6 PfSPZ" | treat == "Placebo") %>%
  mutate(mal.atp.3 = factor(ifelse(mal.atp.3==1, "not protected", "protected"))) %>%
  dplyr::rename(Outcome = "mal.atp.3") %>%
  mutate(mal.vax.1 = ifelse(mal.vax.1==1, "yes", "no")) %>%
  mutate(logCSPIgG = log10(`anti-CSP IgG`+1)) %>%  
   dplyr::rename('Parasitemic at first vax' = "mal.vax.1") %>%
   column_to_rownames(var = "SAMPLEID") %>%
  drop_na(logCSPIgG, `FACS_CD14+_of_live_monocytes`)

Mono_v_CSP_plot_nrow2 <- Mono_v_CSP_plot_nrow2_dat %>%
  ggplot(., aes(x = logCSPIgG, y = `FACS_CD14+_of_live_monocytes`, color = Outcome)) +
  geom_point(size = 1) +
  stat_smooth(method = `lm`, alpha = 0.2, weight = 0.25, size = 0.5) +
  #scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
  ylab("% of parent") +
  xlab("log10(CSP-specific IgG + 1)") +
  scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
  theme_bw(base_family = basefont, base_size = basefontsize) +
  theme(strip.background = element_blank(),
        legend.position="top"
          ) +
  facet_wrap(~treat, nrow = 2, scales = "free_y")

 CSP_pos_by_dose_nrow2_dat <- Mono_v_CSP_plot_nrow2_dat %>%
   dplyr::filter(logCSPIgG >0) %>%
   dplyr::filter(treat == "1.8 x 10^6 PfSPZ" | treat == "Placebo")
 
 CSP_pos_by_dose_nrow2_summarized_dat <- CSP_pos_by_dose_nrow2_dat %>%
   group_by(treat, Outcome) %>%
   summarize(n= n()) %>%
   mutate(n_label = paste0("n = ",n,""))
 
 CSP_pos_by_dose_nrow2_dat <- CSP_pos_by_dose_nrow2_dat %>%
   left_join(., CSP_pos_by_dose_nrow2_summarized_dat,
             by = c("treat", "Outcome")) %>%
   mutate(Outcome = ifelse(Outcome == "protected", "P", "NP")) %>%
   mutate(outcome_n = paste0(Outcome, " (", n, ")"))
 
 CSP_pos_by_dose_plot_nrow2 <- CSP_pos_by_dose_nrow2_dat %>%
   ggplot(., aes(x = Outcome, y = `FACS_CD14+_of_live_monocytes`, fill = Outcome, color = Outcome)) +
   geom_point(position = position_jitterdodge(), size = 1) +
   geom_boxplot(alpha = 0.4, color = "grey30", outlier.shape = NA) +
   ggpubr::stat_compare_means(aes(group = Outcome),
                              label = "p.format", method = "wilcox.test",
                              label.x.npc = "center", vjust = 1,
                              symnum.args = symnum.args, show.legend = FALSE, size = pvalfontsize) +
   scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
   scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
   ylab("% of parent") +
   theme_bw(base_family = basefont, base_size = basefontsize) +
   theme(#axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
          #axis.text.y=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.x=element_blank(),
          strip.background = element_blank(),
          legend.position="top"
          ) + 
  geom_text(aes(label=n_label, y = 0.3), color = "black",size=3,vjust=-0.5, check_overlap = TRUE) +
  facet_wrap(~treat, scale = "free_y", nrow=2) 
 
combinedPlot <- ggarrange(Mono_v_CSP_plot_nrow2, CSP_pos_by_dose_plot_nrow2, ncol = 2, common.legend = TRUE)

Plot relationship between CSP IgG and monocytes (Figure S7B)

Perform correlations between % CD14+ monocytes vs CSP-specific IgG stratified by teratment and outcome (per reviewer request). Shown are coefficients and p values for Pearson’s correlations.

Mono_v_CSP_plot_nrow2_dat %>%
  filter(treat == "1.8 x 10^6 PfSPZ" | treat == "Placebo") %>%
  droplevels() %>%
  drop_na(`FACS_CD14+CD16-_of_live_monocytes`) %>%
  group_by(treat, Outcome) %>%
  dplyr::summarize(r = cor.test(`FACS_CD14+CD16-_of_live_monocytes`,logCSPIgG, method = "pearson")[["estimate"]],
                   p_val = cor.test(`FACS_CD14+CD16-_of_live_monocytes`,logCSPIgG, method = "pearson")[["p.value"]]) %>%
  knitr::kable()
treat Outcome r p_val
Placebo not protected -0.5252408 0.0034364
Placebo protected 0.1462363 0.5754477
1.8 x 10^6 PfSPZ not protected 0.1079542 0.6600058
1.8 x 10^6 PfSPZ protected -0.1767947 0.3875994

Logistic regression to assess three-way interaction

Build model and plot results as flextable

logit_res_all <- glm(Outcome~logCSPIgG*`FACS_CD14+CD16-_of_live_monocytes`*treat, data = Mono_v_CSP_plot_nrow2_dat, family="binomial")
#summary(logit_res_all)

sect_properties <- prop_section(
  page_size = page_size(
    orient = "portrait",
    width = 8, height = 10
  ),
  type = "continuous",
  page_margins = page_mar()
)

my_flex_tab <- tbl_regression(logit_res_all, exponentiate = FALSE) %>%
  add_n() %>%
  as_flex_table()

my_flex_tab

Characteristic

N

log(OR)1

95% CI1

p-value

logCSPIgG

91

-2.2

-4.6, -0.49

0.035

FACS_CD14+CD16-_of_live_monocytes

91

-21

-54, 0.56

0.14

treat

91

Placebo

—

—

1.8 x 10^6 PfSPZ

-1.3

-4.9, 1.6

0.4

logCSPIgG * FACS_CD14+CD16-_of_live_monocytes

91

23

5.7, 47

0.031

logCSPIgG * treat

91

logCSPIgG * 1.8 x 10^6 PfSPZ

2.4

0.34, 5.1

0.042

FACS_CD14+CD16-_of_live_monocytes * treat

91

FACS_CD14+CD16-_of_live_monocytes * 1.8 x 10^6 PfSPZ

23

-2.3, 58

0.14

logCSPIgG * FACS_CD14+CD16-_of_live_monocytes * treat

91

logCSPIgG * FACS_CD14+CD16-_of_live_monocytes * 1.8 x 10^6 PfSPZ

-29

-56, -9.5

0.012

1OR = Odds Ratio, CI = Confidence Interval

Examine and visualize three way interaction using perspective plot.

anova(update(logit_res_all,.~.-treat:logCSPIgG:`FACS_CD14+CD16-_of_live_monocytes`),logit_res_all)
## Analysis of Deviance Table
## 
## Model 1: Outcome ~ logCSPIgG + `FACS_CD14+CD16-_of_live_monocytes` + treat + 
##     logCSPIgG:`FACS_CD14+CD16-_of_live_monocytes` + logCSPIgG:treat + 
##     `FACS_CD14+CD16-_of_live_monocytes`:treat
## Model 2: Outcome ~ logCSPIgG * `FACS_CD14+CD16-_of_live_monocytes` * treat
##   Resid. Df Resid. Dev Df Deviance
## 1        84     118.50            
## 2        83     108.82  1   9.6771
new_placebo_dat <- data.frame("FACS_CD14+CD16-_of_live_monocytes" = rep(seq(0, 0.50, by =0.025),10),
                              logCSPIgG = rep(seq(0,5, by = 0.55),21), check.names = FALSE,
                              treat = rep("Placebo", 210)) %>%
  arrange("FACS_CD14+CD16-_of_live_monocytes", logCSPIgG)
new_placebo_dat$pred <- predict(logit_res_all,newdata=new_placebo_dat)
new_placebo_fit <- new_placebo_dat %>%
  dplyr::rename(monocytes = `FACS_CD14+CD16-_of_live_monocytes`) %>%
  dplyr::select(pred, monocytes, logCSPIgG) %>%
  arrange(monocytes, logCSPIgG)

plot_matrix_placebo <- t(acast(new_placebo_fit, monocytes~logCSPIgG, value.var= "pred"))


new_highdose_dat <- data.frame("FACS_CD14+CD16-_of_live_monocytes" = rep(seq(0, 0.50, by =0.025),10),
                              logCSPIgG = rep(seq(0,5, by = 0.55),21), check.names = FALSE,
                              treat = rep("1.8 x 10^6 PfSPZ", 210)) %>%
  arrange("FACS_CD14+CD16-_of_live_monocytes", logCSPIgG)
new_highdose_dat$pred <- predict(logit_res_all,newdata=new_highdose_dat)
new_highdose_fit <- new_highdose_dat %>%
  dplyr::rename(monocytes = `FACS_CD14+CD16-_of_live_monocytes`) %>%
  dplyr::select(pred, monocytes, logCSPIgG) %>%
  arrange(monocytes, logCSPIgG)

plot_matrix_hd <- t(acast(new_highdose_fit, monocytes~logCSPIgG, value.var= "pred"))

par(mfrow=c(1,2),mai=c(0,0.1,0.2,0)+.02)

my_theta <- 305#310
my_phi <- 20 #20

#placebo
persp(x = as.numeric(rownames(plot_matrix_placebo)), 
  y = as.numeric(colnames(plot_matrix_placebo)), 
  z = plot_matrix_placebo,
  xlab = "CSP-specific IgG",
  ylab = "monocytes",
  zlab = "protection",
  zlim = c(-15,40),
  main = "Placebo",
  #ticktype ='detailed',
  theta = my_theta, 
  phi = 20,
  col = "gray", shade = 0.5)

#1.8x10^6 PfSPZ
persp(x = as.numeric(rownames(plot_matrix_hd)), 
  y = as.numeric(colnames(plot_matrix_hd)), 
  z = plot_matrix_hd,
  xlab = "CSP-specific IgG",
  ylab = "monocytes",
  zlab = "protection",
  zlim = c(-15,40),
  main = "1.8x10^6 PfSPZ",
  #ticktype ='detailed',
  theta = my_theta, 
  phi = my_phi,
  col = "red", shade = 0.5)

Read-in module expression data

#from google drive
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1togaBlNIxiDj16FyXTC-r7Qw0A9cG2az"), path = temp, overwrite = TRUE)
baseline_eset <- readRDS(file = dl$local_path)
MES_eset <- baseline_eset
MES_mat <- Gene2ModuleExpressionScores(MES_eset, module_list = "lowBTMs", summary_stat = median)
myMESdat <- MES_mat %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column(var = "SAMPLEID") %>%
  #dplyr::select(c(SAMPLEID, `INFLAMMATORY/TLR/CHEMOKINES`)) %>%
  right_join(., myDat_monocytes_CSPAb,
             by="SAMPLEID") %>%
  dplyr::filter(treat == "1.8 x 10^6 PfSPZ" | treat == "Placebo") %>%
  mutate(mal.atp.3 = ifelse(mal.atp.3==1, "not protected", "protected")) %>%
  dplyr::rename(Outcome = "mal.atp.3") %>%
  mutate(mal.vax.1 = ifelse(mal.vax.1==1, "yes", "no")) %>%
  mutate(logCSPIgG = log10(`anti-CSP IgG`+1)) %>%
  mutate(CSPIgG_cat = factor(ifelse(logCSPIgG>0, "CSP-specific IgG (+)", "CSP-specific IgG (-)"),
         levels = c("CSP-specific IgG (-)", "CSP-specific IgG (+)"))) %>%
  dplyr::rename('Parasitemic at first vax' = "mal.vax.1") %>%
  column_to_rownames(var = "SAMPLEID")

Baseline expression of innate BTMs in Protected and Not Protected children without CSP-specific IgG at baseline

Comparisons with significant differences p < 0.05 correspond with plots in Figure 7B of revised manuscript. Plots correspond with results in Table S8 revised manuscript.

basefont <- "sans"
basefontsize <- 10
pvalfontsize <- 4

lowBTM_noCSPIgG_Outcome_BoxPlots <- myMESdat %>%
  dplyr::filter(logCSPIgG == 0) %>%
  dplyr::select(treat, Outcome, logCSPIgG, CSPIgG_cat, contains("inflamm") | contains("LPS") | contains("TLR") | contains("innate") | contains("NFkB") | contains("antigen") | contains("lysosome") | contains("activat") & contains("dendritic") | contains("activat") & contains("monocyt") | contains("activat") & contains("DC") | contains("activat") & contains("myeloid") | contains("IFN") | contains("interferon") | contains("viral") | contains("complement")) %>%
  pivot_longer(., cols = c(5:ncol(.)), names_to = "variable", values_to = "values") %>%
  ggplot(., aes(x = treat, y = values, fill = Outcome, color = Outcome)) +
  geom_point(position = position_jitterdodge(), size = 0.25) +
  geom_boxplot(alpha = 0.4, color = "grey30", outlier.shape = NA, weight = 0.5, size = 0.25) +
  ggpubr::stat_compare_means(aes(group = Outcome),
                             label = "p.format", method = "wilcox.test", label.x.npc = "center", vjust = 1,
                             symnum.args = symnum.args, show.legend = FALSE, size = pvalfontsize) +
  scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
  scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
  ylab("module expression score") +
  theme_bw(base_family = basefont, base_size = basefontsize) +
  theme(axis.ticks.x=element_blank(),
        axis.title.x=element_blank(),
        strip.background = element_blank(),
        legend.position="top"
        ) +
  facet_wrap(~variable, scales = "free_y", ncol=4,
             labeller = labeller(variable = label_wrap_gen(width = 25)))

lowBTM_noCSPIgG_Outcome_BoxPlots

Baseline expression of innate BTMs in Protected and Not Protected children with CSP-specific IgG at baseline

Plots correspond with results in Table S8 revised manuscript.

basefont <- "sans"
basefontsize <- 10
pvalfontsize <- 4

lowBTM_withCSPIgG_Outcome_BoxPlots <- myMESdat %>%
  dplyr::filter(logCSPIgG > 0) %>%
  dplyr::select(treat, Outcome, logCSPIgG, CSPIgG_cat, contains("inflamm") | contains("LPS") | contains("TLR") | contains("innate") | contains("NFkB") | contains("antigen") | contains("lysosome") | contains("activat") & contains("dendritic") | contains("activat") & contains("monocyt") | contains("activat") & contains("DC") | contains("activat") & contains("myeloid") | contains("IFN") | contains("interferon") | contains("viral") | contains("complement")) %>%
  pivot_longer(., cols = c(5:ncol(.)), names_to = "variable", values_to = "values") %>%
  ggplot(., aes(x = treat, y = values, fill = Outcome, color = Outcome)) +
  geom_point(position = position_jitterdodge(), size = 0.25) +
  geom_boxplot(alpha = 0.4, color = "grey30", outlier.shape = NA, weight = 0.5, size = 0.25) +
  ggpubr::stat_compare_means(aes(group = Outcome), label = "p.format", method = "wilcox.test", label.x.npc = "center", vjust = 1, symnum.args = symnum.args, show.legend = FALSE, size = pvalfontsize) +
  scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
  scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
  ylab("module expression score") +
  theme_bw(base_family = basefont, base_size = basefontsize) +
  theme(axis.ticks.x=element_blank(),
        axis.title.x=element_blank(),
        strip.background = element_blank(),
        legend.position="top"
        ) +
  facet_wrap(~variable, scales = "free_y", ncol=4,
             labeller = labeller(variable = label_wrap_gen(width = 25)))

lowBTM_withCSPIgG_Outcome_BoxPlots

Plot boxplots between innate low BTM expression in P vs NP stratfied by CSP-specific IgG presence/absence and treatment

Shown is p value for Wilcoxon test between P and NP.

#Figure Wilcoxon test between innate low BTM expression in P vs NP stratefied by CSP-specific IgG presence/absence and treatment
basefont <- "sans"
basefontsize <- 10
pvalfontsize <- 4

addSmallLegend <- function(myPlot, pointSize = 0.25, textSize = 6, spaceLegend = 0.6) {
    myPlot +
        guides(shape = guide_legend(override.aes = list(size = pointSize)),
               color = guide_legend(override.aes = list(size = pointSize))) +
        theme(legend.title = element_text(size = textSize), 
              legend.text  = element_text(size = textSize),
              legend.key.size = unit(spaceLegend, "lines"))
}

lowBTM_innate_splitCSPIgG_Outcome_BoxPlots <- myMESdat %>%
  #dplyr::filter(logCSPIgG > 0) %>%
  dplyr::select(treat, Outcome, logCSPIgG, CSPIgG_cat, contains("inflamm") | contains("LPS") | contains("TLR") | contains("innate") | contains("NFkB") | contains("antigen") | contains("lysosome") | contains("activat") & contains("dendritic") | contains("activat") & contains("monocyt") | contains("activat") & contains("DC") | contains("activat") & contains("myeloid") | contains("IFN") | contains("interferon") | contains("viral") | contains("complement")) %>%
  dplyr::select(!contains("CytokineObs")) %>%
  dplyr::select(!contains("esting")) %>%
  dplyr::select(!contains("ICS")) %>%
  pivot_longer(., cols = c(5:ncol(.)), names_to = "variable", values_to = "values") %>%
  ggplot(., aes(x = treat, y = values, fill = Outcome, color = Outcome)) +
  geom_point(position = position_jitterdodge(), size = 0.25) +
  geom_boxplot(alpha = 0.4, color = "grey30", outlier.shape = NA, weight = 0.5, size = 0.25) +
  ggpubr::stat_compare_means(aes(group = Outcome), label = "p.format", method = "wilcox.test", label.x.npc = "center", vjust = 1, symnum.args = symnum.args, show.legend = FALSE, size = 3) +
  scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
  scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
  ylab("module expression score") +
  theme_bw(base_family = basefont, base_size = basefontsize) +
  theme(#axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        #axis.text.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x=element_blank(),
        strip.background = element_blank(),
        legend.position="top"
        ) +
  facet_wrap(vars(variable, CSPIgG_cat), scales = "free_y", ncol=4,
             labeller = labeller(variable = label_wrap_gen(width = 25)))

lowBTM_innate_splitCSPIgG_Outcome_BoxPlots

Table S8 Wilcoxon test between innate low BTM expression in P vs NP stratified by CSP-specific IgG presence/absence and treatment

Show only comparisons with p<0.05.

lowBTM_innate_splitCSPIgG_Outcome_BoxPlots_data <- ggplot_build(lowBTM_innate_splitCSPIgG_Outcome_BoxPlots)
stat_res <- lowBTM_innate_splitCSPIgG_Outcome_BoxPlots_data$data[[3]]
values_res <- lowBTM_innate_splitCSPIgG_Outcome_BoxPlots_data$data[[2]]
foodat <- myMESdat %>%
  #select only pertinent variables and innate BTMs
  dplyr::select(treat, Outcome, logCSPIgG, CSPIgG_cat, contains("inflamm") | contains("LPS") | contains("TLR") | contains("innate") | contains("NFkB") | contains("antigen") | contains("lysosome") | contains("activat") & contains("dendritic") | contains("activat") & contains("monocyt") | contains("activat") & contains("DC") | contains("activat") & contains("myeloid") | contains("IFN") | contains("interferon") | contains("viral") | contains("complement")) %>%
  dplyr::select(!contains("CytokineObs")) %>%
  dplyr::select(!contains("esting")) %>%
  dplyr::select(!contains("ICS")) %>%
  pivot_longer(., cols = c(5:ncol(.)), names_to = "variable", values_to = "values") %>%
  droplevels()

#Convert to Table S8
res <- foodat %>%
  dplyr::select(-c(logCSPIgG)) %>%
  group_by(variable, CSPIgG_cat, treat) %>% 
  do(w = wilcox.test(values~Outcome, data=., paired=FALSE, exact = TRUE)) %>%
       dplyr::summarise(treat, variable, CSPIgG_cat,
                 Wilcox_stat = w$statistic,
                 Wilcox_p = w$p.value) 
res2 <- foodat %>%
  dplyr::select(-c(logCSPIgG)) %>%
  group_by(variable, CSPIgG_cat, treat, Outcome) %>%
  dplyr::summarise(median = median(values)) %>%
  pivot_wider(., names_from = Outcome, values_from = median)

res_combine <- res2 %>%
  right_join(., res,
            by = c("variable", "CSPIgG_cat", "treat")) %>%
  rename(median_NP = "not protected",
         median_P = "protected") %>%
  arrange(variable)

res_combine %>%
  arrange(Wilcox_p) %>%
  filter(Wilcox_p<0.05) %>%
  knitr::kable()
variable CSPIgG_cat treat median_NP median_P Wilcox_stat Wilcox_p
enriched in antigen presentation (II) (M95.0) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 7.524062 7.7533129 79 0.0014874
lysosome (M209) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.572989 7.3690154 105 0.0025021
complement activation (II) (M112.1) CSP-specific IgG (-) Placebo -1.012838 0.0774692 11 0.0049747
enriched in antigen presentation (III) (M95.1) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.526421 6.1319940 102 0.0053683
complement and other receptors in DCs (M40) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 4.415888 5.0996502 95 0.0073728
viral sensing & immunity; IRF2 targets network (II) (M111.1) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 5.868480 5.4539298 99 0.0105980
enriched in antigen presentation (II) (M95.0) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.870350 7.6310593 97 0.0160241
type I interferon response (M127) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.294684 6.6394908 97 0.0160241
innate antiviral response (M150) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.294684 6.3835160 96 0.0194855
viral sensing & immunity; IRF2 targets network (I) (M111.0) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 5.235057 4.7344488 96 0.0194855
lysosome (M209) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 6.939543 7.2739822 111 0.0278341
MHC-TLR7-TLR8 cluster (M146) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 7.084985 5.6724442 94 0.0282632
antiviral IFN signature (M75) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 5.305919 4.6556607 94 0.0282632
regulation of antigen presentation and immune response (M5.0) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 6.629141 6.8154029 112 0.0299998
inflammasome receptors and signaling (M53) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 6.479248 6.3989044 92 0.0399938
inflammasome receptors and signaling (M53) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 5.983784 6.3304068 118 0.0461582
myeloid, dendritic cell activation via NFkB (I) (M43.0) CSP-specific IgG (+) 1.8 x 10^6 PfSPZ 3.875170 4.2048924 118 0.0461582
regulation of antigen presentation and immune response (M5.0) CSP-specific IgG (-) 1.8 x 10^6 PfSPZ 6.919876 6.8190109 91 0.0471708

Low BTMs that are significantly different (p<0.05) at pre-vaccination baseline between P and NP in infants who without CSP-specific IgG and who would receive 1.8x10^6 PfSPZ

basefont <- "sans"
basefontsize <- 10
pvalfontsize <- 4

#select significant BTMs (p<0.05) in CSP IgG (-) among high-dose PfSPZ group only
res_signif <- res_combine %>%
  filter(Wilcox_p<0.05 & CSPIgG_cat == "CSP-specific IgG (-)" & treat == "1.8 x 10^6 PfSPZ")
lowBTM_innate_significant_CSPneg_highdosePfSPZ_BoxPlots <- myMESdat %>%
  dplyr::filter(logCSPIgG == 0 & treat == "1.8 x 10^6 PfSPZ") %>%
  dplyr::select(treat, Outcome, logCSPIgG, CSPIgG_cat, res_signif$variable) %>%
  mutate(Outcome = ifelse(Outcome == "protected", "P", "NP")) %>%
  pivot_longer(., cols = c(5:ncol(.)), names_to = "variable", values_to = "values") %>%
  ggplot(., aes(x = Outcome, y = values, fill = Outcome, color = Outcome)) +
  geom_point(position = position_jitterdodge(), size = 0.35) +
  geom_boxplot(alpha = 0.4, color = "grey30", outlier.shape = NA, weight = 0.25, size = 0.15) +
  ggpubr::stat_compare_means(aes(group = Outcome), label = "p.format", method = "wilcox.test", label.x.npc = "center", vjust = 1, symnum.args = symnum.args, show.legend = FALSE, size = 3) +
  scale_fill_manual(values = c("#A6CEE3", "#1F78B4")) +
  scale_color_manual(values = c("#A6CEE3", "#1F78B4")) +
  ylab("module expression score") +
  theme_bw(base_family = basefont, base_size = basefontsize) +
  theme(axis.ticks.x=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        strip.background = element_blank(),
        legend.position="top"
        ) +
  facet_wrap(vars(variable), scales = "free_y", ncol=4,
             labeller = labeller(variable = label_wrap_gen(width = 25)))

lowBTM_innate_significant_CSPneg_highdosePfSPZ_BoxPlots